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Abstract. We present results of large-scale three-dimensional simulations of weakly 
magnetized supersonic turbulence at grid resolutions up to 1024"^ cells. Our numerical 
experiments are carried out with the Piecewise Parabolic Method on a Local Stencil and 
assume an isothermal equation of state. The turbulence is driven by a large-scale isotropic 
solenoidal force in a periodic computational domain and fully develops in a few flow crossing 
times. We then evolve the flow for a number of flow crossing times and analyze various 
statistical properties of the saturated turbulent state. We show that the energy transfer rate 
in the inertial range of scales is surprisingly close to a constant, indicating that Kolmogorov's 
phenomenology for incompressible turbulence can be extended to magnetized supersonic flows. 
We also discuss numerical dissipation effects and convergence of different turbulence diagnostics 
as grid resolution refines from 256^ to 1024^ cells. 



1. Introduction 

Supersonic turbulence plays an important role in shaping hierarchical internal substructure of 
molecular clouds (MCs) [H This turbulence is the key ingredient in statistical theories of 
star formation as it controls initial conditions for gravitationally collapsing objects [3 HI E]- 
The nature of molecular cloud turbulence is still poorly understood and no simple conceptual 
theory exists akin to that of Kolmogorov's phenomenology for turbulence in fluids. Two major 
complications are strong compressibility of interstellar gas mediated by radiative energy losses 
and presence of dynamically important magnetic fields. While these two factors bring into play 
additional variables and equations rendering the modeling hardly tractable analytically and 
complicating numerical analysis, the nature of underlying nonlinearity contained in the advection 
terms remains the same. From a mathematician's perspective, the general problem is essentially 
the same as in the Navier-Stokes turbulence, one of the remaining six unsolved Millennium 
Prize Problems for which the Clay Institute of Mathematics is offering a bonus of $1 milhonQ 
Thus hopes to find universality in interstellar turbulence are not entirely unfounded, see e.g. 
[6j, although some researchers believe that breakdown of universality in incompressible isotropic 
magnetohydrodynamic (MHD) turbulence is possible, e.g. due to nonlocality of interactions 



The supersonic regime typical of MC turbulence is extremely hard to achieve in the laboratory 
and the information available from astronomical observations is rather limited [6]. This makes 
numerical simulations, perhaps, the only tool available to explore the statistics of supersonic 
turbulence in detail. High dynamic range simulations can shed light on the energy transfer 
between scales and on the key spatial correlations of the relevant fields in these flows only if they 
provide sufficient scale separation to resolve the dynamics of prevailing nonlinear interactions 
in the inertial range. In addition, low numerical dissipation and high-order numerical methods 
for three-dimensional compressible MHD are needed to achieve an accurate description of the 
inertial range of supersonic turbulence. 

Stable, accurate, divergence-free simulation of magnetized supersonic turbulence is a severe 
test of numerical MHD schemes and has been surprisingly difficult to achieve due to the range of 
flow conditions present. Over the last two years we have developed Piecewise Parabolic Method 
on a Local Stencil (PPML), a new higher order-accurate, low dissipation numerical method which 
requires no additional dissipation or local "fixes" for stable execution [9l 110^ [TT]. PPML is a 
local stencil variant of the popular PPM algorithm [12j for solving the equations of compressible 
ideal magnetohydro dynamics. The principal difference between PPML and PPM is that cell 
interface states are "evolved" rather that reconstructed at every timestep, resulting in a compact 
stencil. Interface states are evolved using Riemann invariants containing all transverse derivative 
information. The scheme is fully multidimensional as the conservation laws are updated in an 
unsplit fashion, including the terms corresponding to the tangential directions in the amplitude 
equations. This helps to avoid numerous well-known pathologies, such as "carbuncles", etc. [13]. 
Divergence-free evolution of the magnetic field is maintained using the higher order-accurate 
constrained transport technique [Hj . The low dissipation and wide spectral bandwidth of this 
method make it an ideal choice for direct simulations of compressible turbulence. 

In Ref. [15] we reported results from a set of three simulations of driven isothermal supersonic 
turbulence at a sonic Mach number of 10 on 512'^ meshes demonstrating the performance of our 
PPML solver on models with different degrees of magnetization from super- Alfvenic through 
trans- Alfvenic regimes. Turbulence statistics derived from these three-dimensional models show 
that the density, velocity, and magnetic energy spectra vary strongly with the strength of 
turbulent magnetic field fluctuations. At the same time, the three cases suggest that a linear 
scaling in the 4/3-law of incompressible MHD turbulence [16^117) is preserved even in the strongly 
compressible regime at Mg = 10, if the proper density weighting is applied to the Elsasser fields 
following the 1/3-rule introduced in Refs. [Ml HSj [23l [19] , see Section |42] for more detail. 

In this paper we present results from PPML simulations of statistically isotropic MHD 
turbulence at a sonic Mach number Mg = 10 and Alfvenic Mach number = 3 on grids 
from 256^ to 1024^ cells. We use this set of simulations to check the (self-)convergence rates for 
various statistical measures with improved grid resolution. We also confirm the universal linear 
scaling for the mixed third-order structure functions of modified Elsasser fields discussed earlier 
in [15]. 



2. Numerical Model and Parameters 

We use PPML to solve numerically the MHD equations for an ideal isothermal gas in a cubic 
domain of size L = 1 with periodic boundary conditions 
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Here p and u are the gas density and velocity, B is the magnetic field, the gas pressure p = c^p 
and the isothermal sound speed = 1, I is the unit tensor. PPML conserves mass, momentum 
and magnetic flux, and keeps V • B = to the machine precision. A fixed large-scale (fc < 2) 
isotropic nonhelical solenoidal force F = pa.— {pa), where a is the acceleration and (...) indicates 
averaging over the entire computational domain, is used to stir the gas keeping Ms close to 10 
during the simulation. The models are initiated with a uniform density po = 1, large-scale 
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velocity field uq = ra, where r is a constant such that (uq) = 8, and a uniform magnetic field 
Bq aligned with the x-coordinate direction. The magnetic field strength is parameterized by 
the ratio of thermal-to- magnetic pressure Pq = pq/Bq'^. In Ref. [15] we carried out three pilot 
simulations with f3o = 0.2, 2, and 20 at 512^. Here we present results for three simulations with 
(3o = 2 and grid resolutions 256^, 512^, and 1024^ cells. The largest of these three experiments 
was carried out using 2048 cores of Ranger supercomputer at TACC. 

Our numerical experiments covered the transition to turbulence and evolution for up to eight 
dynamical times, td = L/2Ms. The magnetic field was not forced and receives energy passively 
through interaction with the velocity field. This random forcing does not generate a mean field, 
but still leads to amplification of the small-scale magnetic energy through compression in shocks 
in combination with flux freezing and, possibly to some extent, via a process known as small-scale 
dynamo, e.g. [20] and references therein. It is implicitly assumed that the effective magnetic 
Prandtl number Pm ~ 1 in our numerical models. This is far from the realistic values for the 
interstellar conditions, where Pm S> 1. The limitations imposed by the available computational 
resources allow us to reach effective Reynolds numbers on the order of 10^ in our largest MHD 
simulations, while the realistic values for molecular clouds can be as large as ~ 10^. 

3. The Saturated Turbulent State and Convergence 

After the initial stir-up phase, our simulations reach a statistical steady state. A snapshot of 
the projected (column) density in such a saturated state is shown in Fig. [TJ The morphology 
of structures visible in this projected density distribution is surprisingly similar to that from 
a non- magnetized simulation at Mg = 6 presented in Fig. 1. of Ref. [15], even though the 
power spectrum index is about —1.7 here (see Sectional and Fig. [3]i) and —2.0 there (see 
Fig. 2a in Ref. [15]) and the morphology of the density distribution in thin slices (not shown) is 
clearly different. This illustrates the loss of information in the projection procedure that makes 
it difficult to measure physical parameters of supersonic turbulent flows from the observed 
projected density distributions alone. 

The properties of this saturated state representing fully developed macroscopically isotropic 
compressible MHD turbulence are the main focus of this paper. While some level of physical 
understanding of small-scale rms field amplification and its saturation exists, the structure of 
the saturated state is poorly known, even in the incompressible case [7]. 

We compared the saturation levels for several integral characteristics of turbulence for models 
with three different resolutions: 256^, 512^, and 1024^. Since in the /3o = 2 case the kinetic energy 
Ek = / \pv?dV dominates over the magnetic energy Em = J \B'^dV by a factor of about 3, 
it is practical to consider Ek and Em separately instead of concentrating on the total energy 
Et = Ek + Em- Figure [2^ shows that the kinetic energy is already converged on 256^ grids, 
in agreement with [21]. In contrast, the saturation level of magnetic energy at /Jq = 2 clearly 
depends on the grid resolution and there are good chances that convergence can be reached 
already at resolution of 2048^ cells. The mean values of Em for t G [3, 7] are 10.0, 12.3, and 13.5 
with the highest corresponding to our finest grid resolution of 1024^. The lack of convergence in 
Em stems from numerical diffusion suppressing small-scale compressions and also the small-scale 
dynamo action, if any, in low-resolution runs. The mean value of Ek for the same averaging 
interval is 39.8, thus the magnetic energy contributes about one quarter of the total energy, 
Em/Et « 0.25, in the 1024^ model. 




Figure 1. Projected (column) density image for a snapshot at t = 6.35td from a 1024^ PPML 
simulation of MHD turbulence at Mg = 10 and Ma = 3 carried out on Ranger (TACC) using 
2048 cores. White-blue-yellow colors correspond to low-intermediate-high projected density 
values. The dynamic range of the image is about 100. 

While the root mean square (rms) sonic Mach number Mg saturates at a level of about 9.2 
independent of the grid resolution, the Alfven Mach number Ma ~ 4.0, 3.1, and 2.8 for grids 
with 256^, 512^, and 1024^ cells, respectively (see Fig. With the same initial uniform 

field Bo, we get turbulent states with slightly different Ma values with a tendency towards less 
super-Alfvenic flows at higher grid resolutions. 

Since we do not include the explicit dissipative terms in the equations we integrate 
numerically, the dissipation in our models is purely numerical and the higher resolution models 
correspond to higher effective Reynolds numberso This can be traced in the non-converging 

We follow here the approach developed in Ref. [22] for nonmagnetized simulations with PPM and assume that 
the nature of dissipation (purely numerical in our case) does not affect the dynamics in the inertial range of scales. 
Resolving the explicit dissipative terms in Eqs. ([2j and ([3j would require a much higher grid resolution to obtain 
an equivalent scale separation. 






Figure 2. Statistics of MHD turbulence at Ms = 10 and /3o — 2, from PPML simulations at grid resolution of 
256^ 512^ and 1024^ cells: (a) mean kinetic, magnetic, and residual energy {gr = [Ek — Em)/ Et) vs. time; (b) 
rms sonic and Alfvenic Mach numbers; (c) mean and rms helicity; (d) mean and rms cross-helicity; (e) rms gas 
density; (f) compensated density power spectra a.t t = 3.5td; (g) compensated velocity power spectra &t t = S.Std; 
(h) compensated magnetic energy power spectra aX t = 3.5td. 



levels of rms helicity gradually growing with higher grid resolution, see Fig. [2]c. Since our driving 
force is non-helical to a very good approximation, the mean helicity stays close to zero over the 
course of a simulation. Figure [2]d shows a rather slow convergence in the levels of the rms 
cross-helicity and finite deviations of the mean cross-helicity Hq = / u • BdV from zero in low 
resolution runs. PPML does not exactly conserve cross-helicity, even though it is one of the ideal 
invariants in MHD. The normalized mean cross-helicity is, however, contained within ±0.8% for 
/So = 2 at 1024^. Even though both the kinetic and the magnetic energy as well as the rms 
Alfven Mach number reach saturation after about 3 dynamical times of evolution in our models 
with /3o = 2, the rms cross-helicity continues to grow slowly for up to 4 or even 5 dynamical 
times meaning that stirring the flow for Std might not be sufficient to reach a statistical steady 
state. 

Another global diagnostic to look at is the rms density. In models with higher resolution the 
flow is better resolved and higher level of density fluctuations is observed. The rms denstity is, 
however, not very sensitive to the grid resolution and changes from 2.2 to 2.3 as the resolution 
jumps from 256^ to 1024'^, see Fig. [2^. 

Finally, to illustrate the convergence of spectral properties of turbulent fluctuations, in Fig.[2]F, 
g, and h we show compensated power spectra of the density, the velocity, and the magnetic field 
for a set of snapshots from simulations with different resolutions taken early in the saturated 
state at t = S.btd- The velocity spectra can be used to access the spectral bandwidth of PPML. 
As one can see from Fig. [2^, the velocity spectrum is resolved up to log^Q k/kmin ~ 1-1 at 256'^ 
and up to log^Q k/kmin ~ 1-4 at 512'^, consistent with the grid refinement by a factor of two. 

4. Exploring Scaling Laws in Supersonic MHD Turbulence 

Here we present time-averaged statistics for our highest resolution (3q = 2 model at 1024^ 
grid cells and compare these with our previous results for non-magnetized (HD) simulations at 
Ms = 6 and grid resolutions of 1024^ and 2048^ cells [l8l [231 [E] . Our MHD data are based 
on 100 full data snapshots equally spaced in time for t G [3,8]frf, i.e. we do the averaging 
over five dynamical times which makes it, perhaps, the best available sample for compressible 
MHD turbulence at this resolution. Our projected density spectra are based on a sample of 750 
projections, including all three projection directions, covering the same time span t £ [3,8]^,^. 
We first begin with the scaling properties of the basic fields, such as the density, velocity, and 
magnetic field and then discuss ideal invariants, such as the total energy and the energy transfer 
rate that may or may not display universal behavior. 

4.I. Basic Fields 

The probability density function (pdf) of the gas density is one of the most important inputs 
for the theories of star formation in molecular clouds as it determines the amount of the dense 
material available for gravitational collapse to form stars [31124^15]. The ability to predict the pdf 
is important in other extreme astrophysical environments, such as, e.g., accretion onto "seed" 
black holes in the first galaxies [25]. It is also known as one of the most robust statistical 
characteristics for compressible isothermal turbulence as its lognormal shape ^26[ [27[ [28l [29] 
gets quickly established in numerical simulations, e.g. [H]. Here we demonstrate the lognormal 
form of the pdf in our /3o = 2 model which gives a perfect fit for a range of ~ 10^ in probability at 
high densities (Fig. [3^) . Compared to the unmagnetized cases at Mg = 6 studied in fL8\ IT5| , the 
lognormal distribution in the MHD case with Ms = 10 is, as expected at higher Mach numbers, 
somewhat wider. If we parametrize the standard deviation cr as a function of the rms sonic Mach 
number, = ln(l + 6^Mf ), we get b = 0.22 and 0.27 as the best fit values for log^o p G [-1, 2.8] 
in the MHD and HD cases, respectively. The mean of the distrubution is also a function of the 
standard deviation, (Inp) = — ^o"^. For more detail on the dependence of the pdf on the value 
of Po see [15] . 



The density power spectrum (Fig. [Sja) scales roughly as k^"^^'^, shallower than the slope of 
— 1 in the HD case at a smaller Mg = 6 [18]. The power spectrum of the logarithm of density 
(Fig. [3}:) scales somewhat steeper, but still shallower than —5/3 found at = 6 [15]. 

The power spectrum of the logarithm of projected gas density roughly follows the scaling /c~^/^, 
which is somewhat steeper than the scaling of the density power spectrum (Fig. EH), but still 
shallower than —2 measured for the HD case at Mg = 6 [K]- The sensitivity of the density 
spectra to the value of /3o is discussed in detail in p^. 

The velocity power spectrum as well as the results of the Helmholtz decomposition of the 
velocity u into the solenoidal (V • = 0) and dilatational (V x Uc = 0) parts, u = + Uc, 
are shown in Fig. [3^. The velocity power index is very close to —3/2 at k/kmin around 30, 
the value proposed for an isotropic incompressible case with turbulent energy equipartition 
\30\ I31j . and the spectrum is dominated by the solenoidal component which follows the same 
scaling. Note that in the HD simulation of Ref. [18] this range of wavenumbers is buried under 
the bottleneck bump which is known to be much weaker, if exists at all, in the MHD models. 
At lower wavenumbers (but still within the inertial range), the velocity spectrum tends to 
get steeper, as should be expected for a strongly supersonic regime. The spectrum for the 
dilatational component of velocity is generally slightly steeper, as was previously reported for 
non-magnetized turbulence [32^ 118]. Figure [3]F compares the fraction of dilatational component 
in the velocity power, Xcik) = P{uc,k)/P(u,k), from our PPML MHD model at /3o = 2 with 
those in our previous PPM = 6 HD runs at resolutions of 1024^ and 2048^ cells [T8l[23| [T5]. 
In both HD runs Xc{k) ~ 1/3 in the inertial range, even though we used different driving forces 
and different parameterizations for PPM diffusion in these simulations. The fact that Xc(^) 
tends to be close to 1/3 in isotropic flows at high Mach numbers can be explained by purely 
geometric considerations. At the same time, in our MHD model Xc(^) ~ 1/4 or even somewhat 
lower. A lower mean level of saturation for Xc{k) in magnetized flows at high Mach numbers can 
be explained by the additional magnetic pressure term in the momentum conservation equation. 
A similar value for Xc{k) was obtained in Ref. [2l] for a model with stronger field (/3o = 0.02) at 
the same grid resolution, cf. [33]. The A;^-compensated power spectra of dilatational velocities in 
Fig. 14 of Ref. ^21j, however, demonstrate an unusual rise of power near the Nyquist frequency 
that we do not observe in our simulations with PPML. Neither do we see such rise in our HD 
models with PPM (Fig. [3^). The origin of this sharp rise in compressions at high wavenumbers 
in simulations with the Athena code [34] is unclear. 



4-2. Ideal Invariants 

Since the total energy Et is a conserved quantity, it is instructive to follow its distribution 
as a function of scale, Exik) = ^ P{p^^^u,k) +P(B,k) , see Fig. [3^. It appears that to a 

good approximation the total energy follows a simple scaling law Exik) ~ fc^/^. It is hard to 
believe, though, that this scaling exponent is universal. We have shown that in non-magnetized 
turbulence the power index for P(p^/^u, k) does depend on the rms Mach number of the flow 
|18j . In MHD models, the scaling of Exik) would also depend on the saturation level of the 
magnetic energy. Nevertheless, it is interesting to see how a combination of the magnetic and 
kinetic energy, each of which do not display a clear extended scaling range (Fig. [3^) add up to 
result in an extended flat stretch of the total energy spectrum. Is this coincidental? Or is this 
a result of self-organization in supersonic super- Alfvenic MHD turbulence? 

A better candidate for universal scaling is the energy transfer rate. Indeed, our HD 
simulations have shown that P{p^/^u,k) follows the Kolmogorov —5/3 scaling at Mg = 6 
|18j . There is a set of exact scaling laws for homogeneous and isotropic incompressible MHD 
turbulence analogous to the 4/5-law of Kolmogorov for ordinary turbulence in neutral fluids 
[35^ [TBI I17j . The MHD laws can be expressed in terms of Elsasser fields, = u it B/y'lvrp 
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Figure 3. Time-average statistics of MHD turbulence at Ms — 10 and /?o = 2, from PPML simulation at grid 
resolution 1024^ cells: (a) the gas density pdf and lognormal fits for logj^Q p G [—1,2.8]; (b) compensated power 
spectrum of the gas density; (c) as b, but for the logarithm density; (d) as b, but for the projected density and 
for the logarithm of projected density; (e) as b, but for the velocity and its solenoidal and dilatational parts; (f) 
ratio of dilatational to total velocity power xc as a function of wavenumber fc; (g) compensated power spectra 
of total, kinetic, and magnetic energy, (h) compensated third-order structure functions Szit) for generalized 
Elsasser fields from 512^ and 1024^ simulations. 



[36], as = (5z'^{i)[5zf{l)f) = -^e^i, where 6z\\{e) = [z(x + ei) - z(x)] • e, d is the space 
dimension, e is a unit vector with arbitrary direction, ei is the displacement vector, the brackets 
denote an ensemble average, and summation over repeated indices is implied. Equivalently, the 
scaling laws can be rewritten in terms of the basic fields (u, B), but in MHD there are no 
separate exact laws for the velocity or the magnetic field alone. It is important not to neglect 
the correlations between the u and B fields or z+ and z~ fields constrained by the invariance 
properties of the equations in turbulent cascade models. In supersonic turbulence, it is equally 
important to properly account for the density-velocity and density-magnetic field correlations. 
In incompressible MHD, the total energy Et and the cross-helicity Hq play the role of ideal 
invariants and, thus, the (total) energy transfer rate 6"^ = (e"^ + e~)/2 and the cross-helicity 
transfer rate e*^ = (e"*" — e~)/2. For vanishing magnetic field B, one recovers the Kolmogorov 
4/5-law, ([5u|[(£)]'^) = —i/bei, where e is the mean rate of kinetic energy transfer [17j . 

Numerical simulations generally confirm these incompressible scalings, although the Reynolds 
numbers are perhaps still too small to reproduce the asymptotic linear behavior with a desired 
precision and the results are sensitive to statistical errors [37 1 132 1 [38]. Often, the absolute value of 
the longitudinal difference is used and still a linear scaling is recovered numerically, while there is 
no rigorous result for the normalization constant in this case. The third-order transverse velocity 
structure functions also show a linear scaling in simulations of nonmagnetized turbulence \32\ [T8] 
and the difference of scaling exponents for longitudinal and transverse structure functions can 
serve as a robust measure of statistical uncertainty of the computed exponents [39]. In MHD 
simulations, the third-order structure functions of the Elsasser fields {^|(52;|j^(^)|^^ were shown to 
have an approximately linear scaling too \37\ HQ] , but again there is no reason for this result to 
universally hold in all situations since the correlations between the fields play an important 
role in nonlinear transfer processes. 

Can the 4/3-law for incompressible MHD turbulence be extended to supersonic regimes in 
molecular clouds? As we discussed in Ref. jl5|, a proper density weighting of the velocity, 
p^/'^u, preserves the approximately linear scaling of the third-order structure functions at high 
Mach numbers in the nonmagnetized case (for a more involved approach to density weighting 
see Ref. [41]). It is straightforward to redefine the Elsasser fields for compressible flows using 
this "1/3-rule", = p^/^{u ± B/^/4^Tp), so that they match the original z^ fields in the 
incompressible limit and reduce to p^^^u in the limit of vanishing B. The new fields can 
have universal scaling properties in homogeneous isotropic turbulent flows with a broad range 
of sonic and Alfvenic Mach numbers. 

We use data from our 1024^ MHD simulation to support the conjecture proposed in Ref. [15] 
and compute the third-order structure functions S^^ and 3 defined above using the fields. 
Since we are interested in the energy transfer through the inertial interval, we compute the sum 
of S^^ and S^^ which determines the energy transfer rate in the incompressible limit. To 
further reduce statistical errors, we combine the transverse and longitudinal structure functions, 
but the absolute value operator was not applied to the field differences. Figure [3h shows the 
compensated scaling for Ss{£) = {S^^ + S^^ + '^'Ja + ^±3)/'^ averaged over 100 flow snapshots 
and the best linear fit Ssii) ~ £i-020±o.ooi f^^, log^/^ g [l.l, 1.8]. For comparison, we overplot 
data from our 512^ simulation discussed in Ref. [12]. Clearly, at 1024^, we get a more extended 
scaling range and confirm the linear scaling found in [15]. Howefer, further work is needed to 
verify this result with high resolution simulations at different levels of magnetization. 

5. Conclusions 

We verified the applicability of the 1/3-rule of hydrodynamic supersonic turbulence to super- 
Alfvenic regimes of MHD turbulence using data from a new simulation with the Piecewise 
Parabolic Method on a Local Stencil on a grid of 1024^ cells. Our results suggest that the 



energy transfer rate is approximately constant accross the inertial interval and that the 4/3-law 
of incompressible MHD can be extended to the supersonic turbulent flows. 
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